O’Hara on GitHub


knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
am_dir <- '~/git-annex/aquamaps_2021'

1 Summary

Read in and examine the newest AquaMaps data.

2 Methods

2.1 Read in and examine species summary file

spp_info <- data.table::fread(file.path(am_dir, 'ver10_2019_speciesoccursum_iucn.csv')) %>%
  janitor::clean_names() %>%
  rename(am_sid = species_id, iucn_sid = iucn_id, comname = f_bname) %>%
  mutate(sciname = tolower(paste(genus, species)))

Number of unique species (by AM species ID): 33518

Number of unique IUCN species (by IUCN species ID): 10519

2.2 Read in and examine species environmental envelope file

env_info <- data.table::fread(file.path(am_dir, 'ver10_2019_hspen.csv')) %>%
  janitor::clean_names() %>%
  rename(am_sid = species_id)

names(env_info)
##  [1] "am_sid"             "speccode"           "fao_areas"         
##  [4] "n_most_lat"         "s_most_lat"         "w_most_long"       
##  [7] "e_most_long"        "depth_yn"           "depth_min"         
## [10] "depth_pref_min"     "depth_pref_max"     "depth_max"         
## [13] "mean_depth"         "pelagic"            "temp_yn"           
## [16] "temp_min"           "temp_pref_min"      "temp_pref_max"     
## [19] "temp_max"           "salinity_yn"        "salinity_min"      
## [22] "salinity_pref_min"  "salinity_pref_max"  "salinity_max"      
## [25] "prim_prod_yn"       "prim_prod_min"      "prim_prod_pref_min"
## [28] "prim_prod_pref_max" "prim_prod_max"      "ice_con_yn"        
## [31] "ice_con_min"        "ice_con_pref_min"   "ice_con_pref_max"  
## [34] "ice_con_max"        "oxy_yn"             "oxy_min"           
## [37] "oxy_pref_min"       "oxy_pref_max"       "oxy_max"           
## [40] "land_dist_yn"       "land_dist_min"      "land_dist_pref_min"
## [43] "land_dist_pref_max" "land_dist_max"      "date_created"      
## [46] "layer"              "rank"               "map_opt"           
## [49] "extn_rule_yn"

Each parameter (depth, temp, salinity, prim_prod, ice_con, oxy, land_dist) has a column for min, pref_min, pref_max, max (and mean for depth). Additionally, all parameters also have columns for _yn (whether that parameter is relevant for that species). For each parameter, let’s plot distributions of the min, max, and preferred min/max for all spp that depend on that parameter.

params <- c('depth', 'temp', 'salinity', 'prim_prod', 'ice_con', 'oxy', 'land_dist')
prm_or <- paste0(params, collapse = '|')
dist_vals <- c('min', 'pref_min', 'mean', 'pref_max', 'max')
env_sum_df <- env_info %>%
  rename(depth_mean = mean_depth) %>%
  select(am_sid, contains(params)) %>%
  gather(param_full, value, -am_sid) %>%
  mutate(param = str_extract(param_full, prm_or),
         dist = str_replace(param_full, paste0('(', prm_or, ')_'), '')) %>%
  select(-param_full) %>%
  spread(dist, value) %>%
  filter(yn == 1) %>%
  select(-yn) %>%
  gather(dist, value, -am_sid, -param) %>%
  mutate(dist = factor(dist, levels = dist_vals),
         value = ifelse(param == 'depth', log10(value + 1), value),
         param = ifelse(param == 'depth', 'log10_depth', param))

write_csv(env_sum_df, file.path(am_dir, 'species_envelope_summary.csv'))

ggplot(env_sum_df) +
  geom_density(aes(x = value, color = dist, ..scaled..), fill = NA) +
  scale_color_viridis_d() +
  facet_wrap(~param, scales = 'free_x') +
  theme_minimal()

env_range_df <- env_sum_df %>%
  mutate(value = ifelse(param == 'log10_depth', 10^value, value)) %>%
  spread(dist, value) %>%
  mutate(mid = ifelse(is.na(mean), (pref_min + pref_max)/2, mean),
         pref_range = pref_max - pref_min,
         max_range  = max - min,
         pref_range = ifelse(param == 'log10_depth', log10(pref_range), pref_range),
         max_range  = ifelse(param == 'log10_depth', log10(max_range), max_range)) %>%
  arrange(pref_range) %>%
  mutate(spp = fct_inorder(am_sid))

ggplot(env_range_df) +
  geom_density(aes(x = pref_range, ..scaled..)) +
  geom_density(aes(x = max_range, ..scaled..), color = 'darkred') +
  facet_wrap(~ param, scales = 'free_x') +
  theme_minimal()

2.3 Read in and examine HCAF file

hcaf_info <- data.table::fread(file.path(am_dir, 'hcaf_v7.csv')) %>%
  janitor::clean_names() %>%
  mutate(across(where(is.numeric), ~ ifelse(.x == -9999, NA, .x)))

names(hcaf_info)
##  [1] "id"                 "csquare_code"       "loiczid"           
##  [4] "n_limit"            "slimit"             "w_limit"           
##  [7] "e_limit"            "center_lat"         "center_long"       
## [10] "cell_area"          "ocean_area"         "p_water"           
## [13] "clim_zone_code"     "fao_area_m"         "fao_area_in"       
## [16] "country_main"       "country_second"     "country_third"     
## [19] "country_sub_main"   "country_sub_second" "country_sub_third" 
## [22] "eez"                "lme"                "lme_border"        
## [25] "meow"               "ocean_basin"        "islands_no"        
## [28] "area0_20"           "area20_40"          "area40_60"         
## [31] "area60_80"          "area80_100"         "area_below100"     
## [34] "elevation_min"      "elevation_max"      "elevation_mean"    
## [37] "elevation_sd"       "depth_min"          "depth_max"         
## [40] "depth_mean"         "depth_sd"           "sst_an_mean"       
## [43] "sbt_an_mean"        "salinity_mean"      "salinity_b_mean"   
## [46] "prim_prod_mean"     "ice_con_ann"        "oxy_mean"          
## [49] "oxy_b_mean"         "land_dist"          "shelf"             
## [52] "slope"              "abyssal"            "tidal_range"       
## [55] "coral"              "estuary"            "seamount"          
## [58] "mpa"
rast_base <- raster(ext = extent(c(-180, 180, -90, 90)), resolution = 0.5, crs = '+init=epsg:4326') %>%
  setValues(1:ncell(.))
var_df <- hcaf_info %>%
  select(loiczid, p_water, eez:land_dist, 
         -starts_with(c('elevation', 'area', 'lme')), -ocean_basin)
vars <- names(var_df)[-1]

for(v in vars) { ### v <- vars[1]
  message('plotting map for ', v)
  var_rast <- subs(rast_base, var_df, by = 'loiczid', which = v)
  plot(var_rast, main = v, col = hcl.colors(n = 100))
}

2.4 Read and examine species cell occurrences

spp_valid <- spp_info %>%
  filter(occur_cells >= 10) %>%
  select(am_sid, sciname)

# csq_loiczid <- hcaf_info %>%
#   select(loiczid, csquare_code) %>%
#   distinct()
# spp_cells <- data.table::fread(file.path(am_dir, 'hcaf_species_native.csv')) %>%
#   janitor::clean_names() %>%
#   left_join(csq_loiczid) %>%
#   select(am_sid = species_id, loiczid, prob = probability)
# 
# write_csv(spp_cells, file.path(am_dir, 'hcaf_species_native_clean.csv'))
spp_cells <- data.table::fread(file.path(am_dir, 'hcaf_species_native_clean.csv')) %>%
  filter(am_sid %in% spp_valid$am_sid)

spp_rich <- spp_cells %>%
  group_by(loiczid) %>%
  summarize(nspp = n_distinct(am_sid), .groups = 'drop')

spp_rich_rast <- subs(rast_base, spp_rich, by = 'loiczid', which = 'nspp')

plot(log10(spp_rich_rast), col = hcl.colors(20))